Intrinsic Activity Develops Along a Sensorimotor-Association Cortical Axis in Youth

Participant Fluctuation Amplitude Data

fluctuations.glasser.pnc <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/PNC/glasserfluctuations_demographics_finalsample.csv") #fMRI + demographics, generated with fitGAMs_fluctuationamplitude_age.R

Glasser Labels

glasser.parcel.labels <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Maps/parcellations/surface/glasser360_regionlist.csv", header = T) #glasser region list

Sensorimotor-Association Axis

S.A.axis.cifti <- read_cifti("/cbica/projects/spatiotemp_dev_plasticity/Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo, vertex-wise axis average ranks
S.A.axis <- as.data.frame(cbind(rank(S.A.axis.cifti$data), names(S.A.axis.cifti$Parcel)))
colnames(S.A.axis) <- c("SA.axis","orig_parcelname")
S.A.axis <- merge(S.A.axis, glasser.parcel.labels, by="orig_parcelname", sort = F)
S.A.axis$SA.axis <- as.numeric(S.A.axis$SA.axis)
rm(S.A.axis.cifti)

Intracortical Myelination Developmental Data

myelin.partialR2.cifti <- read_cifti("/cbica/projects/spatiotemp_dev_plasticity/Myelin/hcpd_n628_myelin_sAge_partial_bayes_r2.pscalar.nii") #T1/T2 ratio - age effect size
myelin.maxdev.cifti <- read_cifti("/cbica/projects/spatiotemp_dev_plasticity/Myelin/hcpd_n628_median_posterior_age_of_max_slope_myelination.pscalar.nii") #age of maximal T1/T2 ratio increase
myelin <- as.data.frame(cbind(myelin.partialR2.cifti$data, myelin.maxdev.cifti$data))
colnames(myelin) <- c("myelin.partialR2","myelin.maxdev")
myelin$label <- glasser.parcel.labels$label
rm(myelin.partialR2.cifti)
rm(myelin.maxdev.cifti)

SNR Mask

SNR.mask <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Maps/parcellations/surface/SNRmask_glasser360.csv")

Spin Test Parcel Rotation Matrix

source("/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/R/rotate.parcellation.R")
source("/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/R/perm.sphere.p.R")
glasser.coords <- read.table("/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/sphere_HCP.txt", header=F) #coordinates of glasser parcel centroids on the freesurfer sphere
perm.id.full <- rotate.parcellation(coord.l = as.matrix(glasser.coords[1:180,]), coord.r = as.matrix(glasser.coords[181:360,]), nrot = 10000) #rotate the glasser parcellation 10,000 times on the freesurfer sphere to generate spatial nulls for spin-based permutation significance testing 
saveRDS(perm.id.full, "/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds")

perm.id.full <- readRDS("/cbica/projects/spatiotemp_dev_plasticity/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins

GAM Results

gam.age.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_age_statistics_glasser.csv") #GAM age smooth statistics, generated with fitGAMs_fluctuationamplitude_age.R

#create df.dev with GAM statistics, S-A axis, plasticity measures, and SNR mask
df.list <- list(gam.age.glasser, S.A.axis, myelin, SNR.mask) #dfs to merge
df.dev <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
df.dev <- df.dev %>% select(label, orig_parcelname, everything()) #order columns
cols = c(3:15)    
df.dev[,cols] = apply(df.dev[,cols], 2, function(x) as.numeric(as.character(x))) #format numerics

#create df.dev.spin formatted for spatial permutation testing
df.dev.spin <- rbind(df.dev[181:360,], df.dev[1:180,]) #format df as left hemisphere -> right hemisphere for spin tests
df.dev.spin$GAM.age.partialR2[df.dev.spin$SNR.mask == 0] <- NA #format data for spin tests by assigning NA to low SNR parcels, treating them like the medial wall
df.dev.spin$minage.decrease[df.dev.spin$SNR.mask == 0] <- NA 

df.dev <- df.dev %>% filter(SNR.mask != 0) #include only high SNR parcels (N=336) in analyses
df.dev$SA.axis.bin <- as.numeric(cut2(df.dev$SA.axis, g=10)) #divide the S-A axis into 10 ranked bins

gam.fitted.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_age_predictedfits_glasser.csv") #GAM predicted fits, generated with fitGAMs_fluctuationamplitude_age.R
gam.fitted.glasser <- inner_join(gam.fitted.glasser, df.dev, by="label", sort = F)

gam.smoothestimates.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_age_smoothestimates_glasser.csv") #GAM smooth estimates, generated with fitGAMs_fluctuationamplitude_age.R
gam.smoothestimates.glasser <- inner_join(gam.smoothestimates.glasser, df.dev, by="label", sort = F)

gam.derivatives.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_age_derivatives_glasser.csv") #GAM true model age smooth derivatives, generated with fitGAMs_fluctuationamplitude_age.R
gam.derivatives.glasser <- left_join(gam.derivatives.glasser, S.A.axis, by="label", sort = F)
gam.derivatives.glasser <- left_join(gam.derivatives.glasser, SNR.mask, by="label", sort = F)
gam.derivatives.glasser <- gam.derivatives.glasser %>% filter(SNR.mask != 0) 

gam.sex.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_agebysex_interaction_glasser.csv") #GAM age spline by sex interaction statistics, generated with fitGAMs_fluctuationamplitude_environment.R

#create df.sex with GAM results, S-A axis, and SNR mask
df.list <- list(gam.sex.glasser, S.A.axis, SNR.mask) #dfs to merge
df.sex <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 

gam.puberty.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_puberty_statistics_glasser.csv") #GAM pubertal stage statistics, generated with fitGAMs_fluctuationamplitude_pubertalstage.R
 
gam.puberty.glasser <- merge(gam.puberty.glasser, SNR.mask, by="label", sort=F)
gam.puberty.glasser <- gam.puberty.glasser %>% filter(SNR.mask != 0) 

gam.env.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_environment_statistics_glasser.csv") #GAM environment statistics, generated with fitGAMs_fluctuationamplitude_environment.R

#create df.env with GAM results, S-A axis, and SNR mask
df.list <- list(gam.env.glasser, S.A.axis, SNR.mask) #dfs to merge
df.env <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 

#create df.env.spin formatted for spatial permutation testing
df.env.spin <- rbind(df.env[181:360,], df.env[1:180,]) #format df as left hemisphere -> right hemisphere for spin tests
df.env.spin$GAM.env.tvalue[df.env.spin$SNR.mask == 0] <- NA #format data for spin tests by assigning NA to low SNR parcels, treating them like the medial wall
df.env.spin$GAM.env.edu.tvalue[df.env.spin$SNR.mask == 0] <- NA #format data for spin tests by assigning NA to low SNR parcels, treating them like the medial wall

df.env <- df.env %>% filter(SNR.mask != 0) #include only high SNR parcels (N=336) in analyses

gam.env.devstage <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_environment_devstage_statistics_glasser.csv") #GAM environment statistics in child, adolescent, and early adult subsamples, generated with fitGAMs_fluctuationamplitude_environment.R

#create df.env.devstage with GAM results, S-A axis, and SNR mask
df.list <- list(gam.env.devstage, S.A.axis, SNR.mask) #dfs to merge
df.env.devstage <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 

#create df.env.devstage.spin formatted for spatial permutation testing
df.env.devstage.spin <- rbind(df.env.devstage[181:360,], df.env.devstage[1:180,]) #format df as left hemisphere -> right hemisphere for spin tests
df.env.devstage.spin$GAM.env.tvalue.child[df.env.devstage.spin$SNR.mask == 0] <- NA #format data for spin tests by assigning NA to low SNR parcels, treating them like the medial wall
df.env.devstage.spin$GAM.env.tvalue.adolescent[df.env.devstage.spin$SNR.mask == 0] <- NA 
df.env.devstage.spin$GAM.env.tvalue.adult[df.env.devstage.spin$SNR.mask == 0] <- NA 

df.env.devstage <- df.env.devstage %>% filter(SNR.mask != 0) #include only high SNR parcels (N=336) in analyses

gam.fitted.env <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_age_predictedfits_byenv.csv") #GAM predicted developmental trajectories for lower and higher factor scores, generated with fitGAMs_fluctuationamplitude_environment.R
gam.fitted.env <- left_join(gam.fitted.env, S.A.axis, by="label", sort = F)
gam.fitted.env <- left_join(gam.fitted.env, SNR.mask, by="label", sort = F)
gam.fitted.env <- gam.fitted.env %>% filter(SNR.mask != 0) 
gam.fitted.env$SA.axis.bin <- as.numeric(cut2(gam.fitted.env$SA.axis, g=10))

Age-dependent changes in intrinsic fMRI activity vary across the cortex

Significant age-related changes in fluctuation amplitude

Number of cortical regions showing significant age effects

sum(p.adjust(df.dev$Anova.age.pvalue, method=c("fdr")) < 0.05) #FDR-corrected p-values from full (+age) versus reduced (-age) models

## [1] 319

Percent of cortical regions showing significant age effects

(sum(p.adjust(df.dev$Anova.age.pvalue, method=c("fdr")) < 0.05))/(nrow(df.dev))*100 

## [1] 94.94048

Across-cortex effect size and direction Figure 1A

ggseg(.data = df.dev, atlas = "glasser", mapping=aes(fill=GAM.age.partialR2), position = c("stacked")) + 
  theme_void() + 
  labs(fill="") +
  paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(-0.07, .03), oob = squish) +
  theme(legend.text = element_text(family = "Arial", color = c("black")))

## merging atlas and data by 'label'

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure1/Brainmap_ageeffect.jpg", device = "pdf", dpi = 500, width = 3 , height = 2)

Age by sex interactions

Number of cortical regions showing significance age*sex interactions

sum(p.adjust(df.sex$GAM.agexsex.pvalue, method=c("fdr")) < 0.05) #FDR-corrected p-values for age by sex interaction. The trajectories of regional age fits did not significantly differ between males and females

## [1] 0

Cortical smooth functions

Figure 1B

gam.smoothestimates.glasser.lh <- gam.smoothestimates.glasser[33401:67200,]
ggplot(gam.smoothestimates.glasser.lh,aes(age,est,group=index,color=GAM.age.partialR2)) + 
  geom_line(size=.3, alpha = .8) + 
  paletteer::scale_color_paletteer_c("pals::ocean.matter", direction = -1, limits = c(-0.07, .03), oob = squish) +
  theme_classic() +
  labs(x = "\nAge", y = "Fluctuation Amplitude (zero centered)\n" ) +
  theme(legend.position = "none") +
  theme(axis.text = element_text(size=6, family = "Arial", color = c("black")), axis.title = element_text(size=6, family = "Arial", color = c("black")), axis.line = element_line(size=.22), axis.ticks = element_line(size=.22)) +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) 

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure1/Regional_agesmooths.jpg", device = "pdf", dpi = 500, width = 2.55 , height = 2.3)

Region-specific developmental trajectories and derivatives

Region V1

V1.dev <- df.dev %>% filter(label == "rh_R_V1")
ggseg(.data = V1.dev, atlas = "glasser", mapping=aes(fill=GAM.age.partialR2, colour=I("#ED8B65"), size=I(.03)), position = c("stacked"), hemisphere = "right", view = "medial" ) + 
  theme_void() + 
  labs(fill="") +
  paletteer::scale_fill_paletteer_c("pals::ocean.matter", na.value="transparent", direction = -1, limits = c(-0.07, .03), oob = squish) +
  theme(legend.position = "none")

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure1/V1_parcel.jpg", device = "pdf", dpi = 500, width = .5 , height = .5)

V1.smooth <- gam.fitted.glasser %>% filter(label == "rh_R_V1")
V1.derivatives <- gam.derivatives.glasser %>% filter(label == "rh_R_V1")

scatterplot <- ggplot(data = fluctuations.glasser.pnc, aes(x = age, y = rh_R_V1)) +
  geom_point(color="#FBDC9D", size = .3) +
  geom_line(data = V1.smooth, aes(x = age, y = fitted), color="#F8B57B", size = 1) +
  geom_ribbon(data = V1.smooth, aes(x = age, y = fitted,  ymin = lower, ymax = upper) ,alpha = .7, linetype = 0, fill="#F8B87D",) +
  labs(x="\nAge", y="Fluctuation Amplitude\n") +
  theme_classic() +
  theme(
  axis.text = element_text(size=6, family = "Arial", color = c("black")),
  axis.title.x=element_text(size=6, family ="Arial", color = c("black")),
  axis.title.y=element_text(size=6, family ="Arial", color = c("black")),
  axis.line = element_line(size=.22), axis.ticks = element_line(size=.22)) +
  scale_x_continuous(breaks=c(8, 12, 16, 20), expand = c(0,.45)) +
  scale_y_continuous(breaks=c(1.0, 2.0), labels = c(sprintf("%.1f", round(1, 2)), sprintf("%.1f", round(2, 2))))

derivplot <- ggplot(data = V1.derivatives) + 
  geom_tile(aes(x = age, y = .1, fill = significant.derivative)) +
  scale_fill_gradient(low = "#EFAF77", high = "#FFE0A1", limits = c(-0.0323371,-.01), na.value = "white") +
  labs(x="\nAge", fill = "Derivative") + 
  theme_classic() +
  theme(legend.position = "none") +
  theme(axis.title.y = element_blank(),
          axis.text = element_blank(),
          axis.line = element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank(),
          legend.text =  element_text(size=6, family = "Arial", color = c("black"))) +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))

allplots <- list(scatterplot,derivplot)
V1.plot <- cowplot::plot_grid(rel_heights = c(16,3.2), plotlist = allplots, align = "v", axis = "lr", ncol = 1)
V1.plot 

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure1/V1_smooth.jpg", device = "pdf", dpi = 500, width = 2.1 , height = 1.85)

Region p24pr

p24pr.dev <- df.dev %>% filter(label == "rh_R_p24pr")
ggseg(.data = p24pr.dev, atlas = "glasser", mapping=aes(fill=GAM.age.partialR2, colour=I("#D54D55"), size=I(.03)), position = c("stacked"), hemisphere = "right", view = "medial" ) + 
  theme_void() + 
  labs(fill="") +
  paletteer::scale_fill_paletteer_c("pals::ocean.matter", na.value="transparent", direction = -1, limits = c(-0.07, .03), oob = squish) +
  theme(legend.position = "none")

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure1/p24pr_parcel.jpg", device = "pdf", dpi = 500, width = .5 , height = .5)

p24pr.smooth <- gam.fitted.glasser %>%filter(label == "rh_R_p24pr")
p24pr.derivatives <- gam.derivatives.glasser %>% filter(label == "rh_R_p24pr")

scatterplot <- ggplot(data = fluctuations.glasser.pnc, aes(x = age, y = rh_R_p24pr)) +
  geom_point(color = "#e06d65", size = .3) +
  geom_line(data = p24pr.smooth, aes(x = age, y = fitted), color="#D54D55", size = 1) +
  geom_ribbon(data = p24pr.smooth, aes(x = age, y = fitted,  ymin = lower, ymax = upper) ,alpha = .7, linetype = 0, fill="#D54D55",) +
  labs(x="\nAge", y="Fluctuation Amplitude\n") +
  theme_classic() +
  theme(
  axis.text = element_text(size=6, family = "Arial", color = c("black")),
  axis.title.x=element_text(size=6, family ="Arial", color = c("black")),
  axis.title.y=element_text(size=6, family ="Arial", color = c("black")),
  axis.line = element_line(size=.22), axis.ticks = element_line(size=.22)) +
  scale_x_continuous(breaks=c(8, 12, 16, 20), expand = c(0,.45)) +
  scale_y_continuous(breaks=c(0.5, 1.0))

derivplot <- ggplot(data = p24pr.derivatives) + 
  geom_tile(aes(x = age, y = .1, fill = significant.derivative)) +
  scale_fill_gradient(high = alpha("#DF5E54",0.7), low = "#D54D55", na.value = "white", limits = c(min(p24pr.derivatives$significant.derivative),-.001)) +
  labs(x="\nAge", fill = "Derivative") + 
  theme_classic() +
  theme(legend.position = "none") +
  theme(axis.title.y = element_blank(),
          axis.text = element_blank(),
          axis.line = element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank(),
          legend.text =  element_text(size=6, family = "Arial", color = c("black"))) +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))

allplots <- list(scatterplot,derivplot)
p24pr.plot <- cowplot::plot_grid(rel_heights = c(16,3.2), plotlist = allplots, align = "v", axis = "lr", ncol = 1)
p24pr.plot 

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure1/p24pr_smooth.jpg", device = "pdf", dpi = 500, width = 2.1 , height = 1.85)  

Region IFSa

IFSa.dev <- df.dev %>% filter(label == "lh_L_IFSa")
ggseg(.data = IFSa.dev, atlas = "glasser", mapping=aes(fill=GAM.age.partialR2, colour=I("#381043"), size=I(.03)), position = c("stacked"), hemisphere = "left", view = "lateral" ) + 
  theme_void() + 
  labs(fill="") +
  paletteer::scale_fill_paletteer_c("pals::ocean.matter", na.value="transparent", direction = -1, limits = c(-0.07, .03), oob = squish) +
  theme(legend.position = "none")

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure1/IFSa_parcel.jpg", device = "pdf", dpi = 500, width = .5 , height = .5)

IFSa.smooth <- gam.fitted.glasser %>%filter(label == "lh_L_IFSa")
IFSa.derivatives <- gam.derivatives.glasser %>% filter(label == "lh_L_IFSa")

scatterplot <- ggplot(data = fluctuations.glasser.pnc, aes(x = age, y = lh_L_IFSa)) +
  geom_point(color="#893774", size = .3) +
  geom_line(data = IFSa.smooth, aes(x = age, y = fitted), color="#381043", size = 1) +
  geom_ribbon(data = IFSa.smooth, aes(x = age, y = fitted,  ymin = lower, ymax = upper) ,alpha = .7, linetype = 0, fill="#381043",) +
  labs(x="\nAge", y="Fluctuation Amplitude\n") +
  theme_classic() +
  theme(
  axis.text = element_text(size=6, family = "Arial", color = c("black")),
  axis.title.x=element_text(size=6, family ="Arial", color = c("black")),
  axis.title.y=element_text(size=6, family ="Arial", color = c("black")),
  axis.line = element_line(size=.22), axis.ticks = element_line(size=.22)) +
  scale_x_continuous(breaks=c(8, 12, 16, 20), expand = c(0,.45)) +
  scale_y_continuous(breaks=c(1.0, 2.0), labels = c(sprintf("%.1f", round(1, 2)), sprintf("%.1f", round(2, 2))))

derivplot <- ggplot(data = IFSa.derivatives) + 
  geom_tile(aes(x = age, y = .1, fill = significant.derivative)) +
  scale_fill_gradient2(high = "#6E195E", low = "#cf93c4", midpoint = 0, mid = "white", na.value = "white") +
  labs(x="\nAge", fill = "Derivative") + 
  theme_classic() +
  theme(legend.position = "none") +
  theme(axis.title.y = element_blank(),
          axis.text = element_blank(),
          axis.line = element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank(),
          legend.text =  element_text(size=6, family = "Arial", color = c("black"))) +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))

allplots <- list(scatterplot,derivplot)
IFSa.plot <- cowplot::plot_grid(rel_heights = c(16,3.2), plotlist = allplots, align = "v", axis = "lr", ncol = 1)
IFSa.plot 

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure1/IFSa_smooth.jpg", device = "pdf", dpi = 500, width = 2.1 , height = 1.85)  

Pubertal stage effects

Association between linear effects of pubertal stage (ordinal variable) and age effects

cor.test(gam.puberty.glasser$GAM.puberty.linear.tvalue, df.dev$GAM.age.partialR2)

## 
##  Pearson's product-moment correlation
## 
## data:  gam.puberty.glasser$GAM.puberty.linear.tvalue and df.dev$GAM.age.partialR2
## t = 8.3455, df = 334, p-value = 1.9e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3227336 0.5001501
## sample estimates:
##       cor 
## 0.4153843

Number of regions where pubertal stage explains significant variance above and beyond age

sum(p.adjust(gam.puberty.glasser$Anova.puberty.pvalue, method=c("fdr")) < 0.05)

## [1] 0

Number of regions with significant age effects when controlling for pubertal stage

gam.age.puberty.glasser <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/fluctuationamplitude_age_pubertycontrolled_statistics_glasser.csv") 
gam.age.puberty.glasser <- merge(gam.age.puberty.glasser, SNR.mask, by="label", sort=F)
gam.age.puberty.glasser <- gam.age.puberty.glasser %>% filter(SNR.mask != 0) 

sum(p.adjust(gam.age.puberty.glasser$Anova.age.pvalue, method=c("fdr")) < 0.05)

## [1] 283

Association between age effects when controlling and not controlling for pubertal stage

cor.test(df.dev$GAM.age.partialR2, gam.age.puberty.glasser$GAM.age.partialR2)

## 
##  Pearson's product-moment correlation
## 
## data:  df.dev$GAM.age.partialR2 and gam.age.puberty.glasser$GAM.age.partialR2
## t = 35.165, df = 334, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8621815 0.9081022
## sample estimates:
##       cor 
## 0.8873217

The development of intrinsic fMRI activity is linked to the maturation of a key regulator of plasticity

Cortical distribution of age effects

Fluctuation amplitude age effects map Figure 2A

ggseg(.data = df.dev, atlas = "glasser", mapping=aes(fill=GAM.age.partialR2), position = c("stacked"), hemisphere = c("right")) + 
  theme_void() + 
  labs(fill="") +
  paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(-0.07, .03), oob = squish) +
  theme(legend.text = element_blank())

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure2/Brainmap_function.jpg", device = "pdf", dpi = 500, width = 2.4, height = 1.7)

T1/T2 ratio age effects map Figure 2A

ggseg(.data = df.dev, atlas = "glasser", mapping=aes(fill=myelin.partialR2), position = c("stacked"), hemisphere = "right") + 
  theme_void() + 
  labs(fill="") +
  paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = 1, limits = c(0.0,0.3), oob = squish) +
  theme(legend.text = element_blank())

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure2/Brainmap_myelin.jpg", device = "pdf", dpi = 500, width = 2.4, height = 1.7)

Correlation (rho + p.spin) between fluctuation amplitude age effects and T1/T2 ratio age effects

cor.test(df.dev$GAM.age.partialR2, df.dev$myelin.partialR2, method=c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  df.dev$GAM.age.partialR2 and df.dev$myelin.partialR2
## S = 10557746, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.6699693

perm.sphere.p(df.dev.spin$GAM.age.partialR2, df.dev.spin$myelin.partialR2, perm.id.full,corr.type='spearman')

## [1] 0.00045

Correlation plot Figure 2B

ggplot(df.dev, aes(x = myelin.partialR2, y = GAM.age.partialR2, fill = GAM.age.partialR2)) + 
geom_point(aes(color = GAM.age.partialR2), shape = 21, size = 1.5) +
paletteer::scale_color_paletteer_c("pals::ocean.matter", direction = -1, limits = c(-0.07, .03), oob = squish) +
paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(-0.07, .03), oob = squish) +
labs(x="\nMyelin Age Effect", y="Fluctuation Amplitude Age Effect\n") +
geom_smooth(method = 'lm', se = TRUE, fill = alpha(c("gray70"),.7), col = "black", size = .75) +
theme_classic() + 
theme(legend.position = "none") +
theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .22), axis.ticks = element_line(size = .22))

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure2/Ageeffects_correlationplot.jpg", device = "pdf", dpi = 500, width = 2.68 , height = 2.175)

Temporal coordination

Fluctuation amplitude age of decrease onset Figure 2C

ggseg(.data = df.dev, atlas = "glasser", mapping = aes(fill = minage.decrease), position = c("stacked"), hemisphere = "right") + 
  theme_void() + 
  labs(fill = "") +
  paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(10, 17.5), oob = squish) +
  theme(legend.text = element_blank())

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure2/Brainmap_function_agedecrease.jpg", device = "pdf", dpi = 500, width = 2.4, height = 1.7)

T1/T2 ratio age of maximal increase Figure 2C

ggseg(.data = df.dev, atlas = "glasser", mapping = aes(fill = myelin.maxdev), position = c("stacked"), hemisphere = "right") + 
  theme_void() + 
  labs(fill = "") +
  paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(10,17.5), oob = squish) +
  theme(legend.text = element_blank())

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure2/Brainmap_myelin_ageincrease.jpg", device = "pdf", dpi = 500, width = 2.4, height = 1.7)

Correlation (rho + p.spin) between fluctuation amplitude age of decrease onset and T1/T2 ratio max age of increase

#All regions
cor.test(df.dev$minage.decrease, df.dev$myelin.maxdev, method=c("spearman"), exact=F)

## 
##  Spearman's rank correlation rho
## 
## data:  df.dev$minage.decrease and df.dev$myelin.maxdev
## S = 1677077, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6382729

perm.sphere.p(df.dev.spin$minage.decrease, df.dev.spin$myelin.maxdev, perm.id.full,corr.type='spearman')

## [1] 0.00125

#Age of decrease onset > 8.5
late.decrease <- df.dev %>% filter(minage.decrease > 8.5)
cor.test(late.decrease$minage.decrease, late.decrease$myelin.maxdev, method=c("spearman"), exact=F)

## 
##  Spearman's rank correlation rho
## 
## data:  late.decrease$minage.decrease and late.decrease$myelin.maxdev
## S = 349670, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6402458

late.decrease.spin <- df.dev.spin
late.decrease.spin$minage.decrease[late.decrease.spin$minage.decrease < 8.5] <- NA
perm.sphere.p(late.decrease.spin$minage.decrease, late.decrease.spin$myelin.maxdev, perm.id.full,corr.type='spearman')

## [1] 0.01565

Correlation plot Figure 2D

ggplot(late.decrease, aes(x = myelin.maxdev, y = minage.decrease, fill = minage.decrease)) + 
geom_point(aes(color = minage.decrease), shape = 21, size = 1.2) +
paletteer::scale_fill_paletteer_c("pals::ocean.matter", direction = -1, limits = c(10, 17.5), oob = squish) +
paletteer::scale_color_paletteer_c("pals::ocean.matter", direction = -1, limits = c(10, 17.5), oob = squish) +
labs(x="\nMyelin Age Effect", y="Fluctuation Amplitude Age Effect\n") +
geom_smooth(method = 'lm', se = TRUE, fill = alpha(c("gray70"),.7), col = "black", size = .75) +
theme_classic() + 
theme(legend.position = "none") +
theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .22), axis.ticks = element_line(size = .22)) +
scale_y_continuous(breaks=c(10, 12, 14, 16, 18, 20))

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure2/Temporaleffects_correlationplot.jpg", device = "pdf", dpi = 500, width = 2.56 , height = 2.175)

Spatiotemporal variability in developmental trajectories is captured by the sensorimotor-association axis

Sensorimotor-association axis Figure 3A

ggseg(.data = df.dev, atlas = "glasser", mapping = aes(fill = SA.axis), position = c("stacked"), hemisphere = "right") + 
  theme_void() + 
  labs(fill = "") +
  scale_fill_gradient2(low = "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 180) +
  theme(legend.text = element_blank())

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure3/SAaxis.jpg", device = "pdf", dpi = 500, width = 3.65, height = 2.95)

Age effects along the S-A Axis

Correlation (rho + p.spin) between fluctuation amplitude age effects and sensorimotor-association axis rank

cor.test(df.dev$GAM.age.partialR2, df.dev$SA.axis, method=c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  df.dev$GAM.age.partialR2 and df.dev$SA.axis
## S = 2934378, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5358554

perm.sphere.p(df.dev.spin$GAM.age.partialR2, df.dev.spin$SA.axis, perm.id.full,corr.type='spearman')

## [1] 0.00215

Temporal development along the S-A Axis

Correlation (rho + p.spin) between fluctuation amplitude age of decrease onset and sensorimotor-association axis rank

cor.test(df.dev$minage.decrease, df.dev$SA.axis, method=c("spearman"), exact=F)

## 
##  Spearman's rank correlation rho
## 
## data:  df.dev$minage.decrease and df.dev$SA.axis
## S = 2477911, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4655417

perm.sphere.p(df.dev.spin$minage.decrease, df.dev.spin$SA.axis, perm.id.full,corr.type='spearman')

## [1] 0.0388

cor.test(late.decrease$minage.decrease, late.decrease$SA.axis, method=c("spearman"), exact=F)

## 
##  Spearman's rank correlation rho
## 
## data:  late.decrease$minage.decrease and late.decrease$SA.axis
## S = 314027, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6769174

perm.sphere.p(late.decrease.spin$minage.decrease, late.decrease.spin$SA.axis, perm.id.full,corr.type='spearman')

## [1] 0.0011

Sensorimotor-association axis binned GAM smooths Figure 3C

gam.smoothestimates.glasser$age <- round(gam.smoothestimates.glasser$age, 3)

meansmooth.bins <- gam.smoothestimates.glasser %>% group_by(age, SA.axis.bin) %>% do(est.mean = mean(.$est)) %>% unnest(cols = c(est.mean)) #calculate the average smooth estimate at each age within each of the 10 S-A axis bins

ggplot(data = meansmooth.bins, aes(x = age, y = est.mean, group = SA.axis.bin, color = as.factor(SA.axis.bin))) +
  geom_line(size = 1.5) +
  scale_color_manual(values = c("#FFC228", "#FFCA4E", "#FFD16A", "#FFDA89", "#FFE6B0", "#D7BCDA", "#BE93C4", "#9859A4", "#863E95", "#6F1382")) +
  theme_classic() +
  labs(x = "\nAge", y = "Fluctuation Amplitude (zero centered)\n", color = "S-A axis deciles") +
  theme(legend.position = c(.585, 0.94), legend.direction = "horizontal", legend.text = element_text(size = 6), legend.key.size = unit(.3, "cm"), legend.title = element_text(size = 6, family = "Arial", color = c("black"))) +
  theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size=.22), axis.ticks = element_line(size = .22)) + scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure3/SAaxis_bin_smooths.jpg", device = "pdf", dpi = 500, width = 3.34 , height = 2.15)

Principal developmental gradient

Principal component analysis to identify the spatial axis that captures the most variance in regional fluctuation amplitude developmental trajectories

gam.smoothestimates.glasser <- gam.smoothestimates.glasser %>% select(label, age, est)
gam.smoothestimates.long <- gam.smoothestimates.glasser %>% pivot_wider(names_from="label",values_from="est")
gam.smoothestimates.long <- gam.smoothestimates.long %>% select(-age) #ages by regions matrix for PCA

smooths.pca <- prcomp(gam.smoothestimates.long) #PCA

Variance explained by PC1

summary(smooths.pca)$importance[2,1]

## [1] 0.87026

Developmental PC1 correlation with the sensorimotor-association axis

PC1 <- as.data.frame(smooths.pca$rotation[,1]) #first principal component loadings
PC1$label <- row.names(PC1) 
colnames(PC1) <- c("PC1","label")
PC1$ranked <- rank(PC1$PC1)

cor.test(PC1$PC1, df.dev$SA.axis, method=c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  PC1$PC1 and df.dev$SA.axis
## S = 1896048, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.700093

df.dev.spin <- left_join(df.dev.spin, PC1, by="label")
perm.sphere.p(df.dev.spin$PC1, df.dev.spin$SA.axis, perm.id.full,corr.type='spearman')

## [1] 0

Correlation plot Figure 3B

PC1 <- left_join(PC1, S.A.axis, by="label")
ggplot(PC1, aes(x = SA.axis, y = PC1, fill = SA.axis)) + 
geom_point(aes(color = SA.axis), shape = 21, size = 1.2) +
scale_fill_gradient2(low = "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 180) +
scale_fill_gradient2(low = "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, midpoint = 180) +
labs(x="\nSensorimtor-Association Axis Rank", y="Principal Developmental Axis Loading\n") +
geom_smooth(method = 'lm', se = TRUE, fill = alpha(c("gray70"),.7), col = "black", size = .25) +
theme_classic() + 
theme(legend.position = "none") +
theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .22), axis.ticks = element_line(size = .22))

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure3/AxisPC1_correlationplot.jpg", device = "pdf", dpi = 500, width = 3.34 , height = 2.15)

Principal developmental axis brain map Figure 3A

ggseg(.data = PC1, atlas = "glasser", mapping=aes(fill = PC1), position = c("stacked"), hemisphere = "right") + 
theme_void() + 
labs(fill = "") +
scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.05, limits = c(-0.11,0.008), oob = squish) + theme(legend.text = element_blank())

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure3/PC1.jpg", device = "pdf", dpi = 500, width = 3.65, height = 2.95)

Developmental PC1 correlation with anatomical, functional, and evolutionary cortical hierarchies

maps <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/Maps/S-A_ArchetypalAxis/Glasser360_MMP_WholeBrain/brainmaps_glasser.csv")
maps$label <- glasser.parcel.labels$label
maps <- merge(PC1, maps, by = "label")

Anatomical hierarchy

cor.test(maps$PC1, maps$T1T2ratio, method = c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  maps$PC1 and maps$T1T2ratio
## S = 10161350, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.6072694

cocor.dep.groups.overlap(
  r.jk = (cor.test(PC1$PC1, df.dev$SA.axis, method=c("spearman"))$estimate), 
  r.jh = (abs(cor.test(maps$PC1, maps$T1T2ratio, method=c("spearman"))$estimate)), 
  r.kh = (abs(cor.test(maps$T1T2ratio, maps$SA.axis, method=c("spearman"))$estimate)), 
  n = 336, alternative = "two.sided", test = "hittner2003")

## 
##   Results of a comparison of two overlapping correlations based on dependent groups
## 
## Comparison between r.jk = 0.7001 and r.jh = 0.6073
## Difference: r.jk - r.jh = 0.0928
## Related correlation: r.kh = 0.7785
## Group size: n = 336
## Null hypothesis: r.jk is equal to r.jh
## Alternative hypothesis: r.jk is not equal to r.jh (two-sided)
## Alpha: 0.05
## 
## hittner2003: Hittner, May, and Silver's (2003) modification of Dunn and Clark's z (1969) using a backtransformed average Fisher's (1921) Z procedure
##   z = 3.5208, p-value = 0.0004
##   Null hypothesis rejected

Functional hierarchy

cor.test(maps$PC1, maps$G1.fMRI, method = c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  maps$PC1 and maps$G1.fMRI
## S = 2560212, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.595039

cocor.dep.groups.overlap(
  r.jk = (cor.test(PC1$PC1, df.dev$SA.axis, method=c("spearman"))$estimate), 
  r.jh = (abs(cor.test(maps$PC1, maps$G1.fMRI, method=c("spearman"))$estimate)), 
  r.kh = (abs(cor.test(maps$G1.fMRI, maps$SA.axis, method=c("spearman"))$estimate)), 
  n = 336, alternative = "two.sided", test = "hittner2003")

## 
##   Results of a comparison of two overlapping correlations based on dependent groups
## 
## Comparison between r.jk = 0.7001 and r.jh = 0.595
## Difference: r.jk - r.jh = 0.1051
## Related correlation: r.kh = 0.8664
## Group size: n = 336
## Null hypothesis: r.jk is equal to r.jh
## Alternative hypothesis: r.jk is not equal to r.jh (two-sided)
## Alpha: 0.05
## 
## hittner2003: Hittner, May, and Silver's (2003) modification of Dunn and Clark's z (1969) using a backtransformed average Fisher's (1921) Z procedure
##   z = 5.0048, p-value = 0.0000
##   Null hypothesis rejected

Evolutionary hierarchy

cor.test(maps$PC1, maps$Evolution.Expansion, method = c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  maps$PC1 and maps$Evolution.Expansion
## S = 4298202, p-value = 2.362e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3201328

cocor.dep.groups.overlap(
  r.jk = (cor.test(PC1$PC1, df.dev$SA.axis, method=c("spearman"))$estimate), 
  r.jh = (abs(cor.test(maps$PC1, maps$Evolution.Expansion, method=c("spearman"))$estimate)), 
  r.kh = (abs(cor.test(maps$Evolution.Expansion, maps$SA.axis, method=c("spearman"))$estimate)), 
  n = 336, alternative = "two.sided", test = "hittner2003")

## 
##   Results of a comparison of two overlapping correlations based on dependent groups
## 
## Comparison between r.jk = 0.7001 and r.jh = 0.3201
## Difference: r.jk - r.jh = 0.38
## Related correlation: r.kh = 0.6026
## Group size: n = 336
## Null hypothesis: r.jk is equal to r.jh
## Alternative hypothesis: r.jk is not equal to r.jh (two-sided)
## Alpha: 0.05
## 
## hittner2003: Hittner, May, and Silver's (2003) modification of Dunn and Clark's z (1969) using a backtransformed average Fisher's (1921) Z procedure
##   z = 9.6478, p-value = 0.0000
##   Null hypothesis rejected

Development is hierarchical through adolescence

Age-specific fluctuation amplitude derivatives

Regional derivative by age plot

gam.derivatives.glasser.lh <- gam.derivatives.glasser[33401:67200,]
ggplot(gam.derivatives.glasser.lh,aes(age,derivative,group=index,color=SA.axis)) + 
  geom_line(size=1, alpha = .7) + 
  scale_color_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, midpoint = 180) +
  theme_classic() +
  theme(legend.position = "none") +
  labs(x = "\nAge", y = "Derivative\n", color = "Sensorimotor-Association Axis Bin") +
  theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .22), axis.ticks = element_line(size = .22)) +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))

gam.derivatives.glasser$age <- round(gam.derivatives.glasser$age, 5)

Age 10 derivative map

age10.derivs <- gam.derivatives.glasser %>% filter(age == 10.03015)

ggseg(.data = age10.derivs, atlas = "glasser", mapping=aes(fill=derivative), position = c("stacked"), hemisphere = "right") + 
  theme_void() + 
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish)

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure4/Age10_derivatives.jpg", device = "pdf", dpi = 500, width = 2.95 , height = 1.95)

Age 10 derivative-axis correlation plot

ggplot(age10.derivs, aes(x = SA.axis, y = derivative, fill = derivative)) + 
geom_point(aes(color = derivative), shape = 21, size = .1) +
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish) +
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish) +
labs(x="\nSensorimtor-Association Axis Rank", y="Derivative\n") +
geom_smooth(method = 'lm', se = TRUE, fill = alpha(c("gray70"),.7), col = "black", size = .35) +
theme_classic() + 
theme(legend.position = "none") +
theme(axis.text = element_blank(), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .15), axis.ticks = element_line(size = .15)) +
scale_y_continuous(limits = c(-.05, .055), oob = squish)

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure4/Age10_correlationplot.jpg", device = "pdf", dpi = 500, width = 1.45 , height = 1.25)

Age 15 derivative map

age15.derivs <- gam.derivatives.glasser %>% filter(age == 15.02429)

ggseg(.data = age15.derivs, atlas = "glasser", mapping=aes(fill=derivative), position = c("stacked"), hemisphere = "right") + 
  theme_void() + 
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish) 

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure4/Age15_derivatives.jpg", device = "pdf", dpi = 500, width = 2.95 , height = 1.95)

Age 15 derivative-axis correlation plot

ggplot(age15.derivs, aes(x = SA.axis, y = derivative, fill = derivative)) + 
geom_point(aes(color = derivative), shape = 21, size = .1) +
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish) +
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish) +
labs(x="\nSensorimtor-Association Axis Rank", y="Derivative\n") +
geom_smooth(method = 'lm', se = TRUE, fill = alpha(c("gray70"),.7), col = "black", size = .35) +
theme_classic() + 
theme(legend.position = "none") +
theme(axis.text = element_blank(), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .15), axis.ticks = element_line(size = .15)) +
scale_y_continuous(limits = c(-.05, .055), oob = squish)

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure4/Age15_correlationplot.jpg", device = "pdf", dpi = 500, width = 1.45 , height = 1.25)

Age 20 derivative map

age20.derivs <- gam.derivatives.glasser %>% filter(age == 20.01843)

ggseg(.data = age20.derivs, atlas = "glasser", mapping=aes(fill=derivative), position = c("stacked"), hemisphere = "right") + 
  theme_void() + 
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.012, limits = c(-0.025,0.005), oob = squish)

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure4/Age20_derivatives.jpg", device = "pdf", dpi = 500, width = 2.95 , height = 1.95)

Age 20 derivative-axis correlation plot

ggplot(age20.derivs, aes(x = SA.axis, y = derivative, fill = derivative)) + 
geom_point(aes(color = derivative), shape = 21, size = .1) +
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish) +
  scale_fill_gradient2(low= "goldenrod1", mid = "white", high = "#6f1282", guide = "colourbar", aesthetics = "color", name = NULL, midpoint = -0.012, limits = c(-0.025, 0.005), oob = squish) +
labs(x="\nSensorimtor-Association Axis Rank", y="Derivative\n") +
geom_smooth(method = 'lm', se = TRUE, fill = alpha(c("gray70"),.7), col = "black", size = .35) +
theme_classic() + 
theme(legend.position = "none") +
theme(axis.text = element_blank(), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .15), axis.ticks = element_line(size = .15)) +
scale_y_continuous(limits = c(-.05, .055), oob = squish)

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure4/Age20_correlationplot.jpg", device = "pdf", dpi = 500, width = 1.45 , height = 1.25)

Regional developmental change across the sensorimotor-association axis Figure 4A

derivative.colorbar <- c("#ff7e38", "#f29049", "#ffa25e", "#faab73", "#fabb73", "#ffcd94", "#ffd894", "#ffe6ab", "#ffe6ab", "#fff5de", "#ffffff", "#f79a97", "#f58682", "#EE6E69", "#C95665", "#BB4166", "#B72863", "#941348", "#7d0f65",  "#660c5f", "#4f0e4a")

ggplot(data=gam.derivatives.glasser,aes(x = age, y = SA.axis, group = index)) +
  geom_line(aes(color=derivative), size=1) +
  scale_color_gradientn(colours = derivative.colorbar, limits=c(-0.028, 0.028), oob = squish) +
  theme_classic() +
  ylab("Sensorimotor-Association Axis Rank\n") +
  xlab("\nAge") +
  labs(color="") +
   theme(axis.text = element_text(size = 7, family = "Arial", color = c("black")), axis.title = element_text(size = 7, family = "Arial", color = c("black")), axis.line = element_line(size =.22), axis.ticks = element_line(size = .22), legend.text = element_text(size = 6, family = "Arial", color = c("black"))) +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45))   

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure4/Derivative_rank_age_plot.pdf", device = "pdf", dpi = 500, width = 3.7 , height = 4.9)

Age-specific correlations between developmental change and the S-A Axis

Posterior derivative - sensorimotor-association axis correlations, N=10,000

deriv.SAaxis.posteriorcorrs <- read.table("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/SAaxis_posteriorderivative_correlation_byage_glasser.csv", header = F, sep = ",") #get the regional derivative-regional axis rank correlation at each age for all draws from the posterior distribution
colnames(deriv.SAaxis.posteriorcorrs) <- sprintf("draw%s",seq(from = 1, to = ncol(deriv.SAaxis.posteriorcorrs)))
deriv.SAaxis.posteriorcorrs <- cbind(deriv.SAaxis.posteriorcorrs, (gam.derivatives.glasser %>% group_by(age) %>% group_keys()))

Age of maximal sensorimotor-association axis correlation: posterior median value + 95% credible interval

age.max.corr <- deriv.SAaxis.posteriorcorrs %>% #find the age at which the correlation is largest for each draw
    summarise(across(contains("draw"),
                     .fns = function(x){
                       round(deriv.SAaxis.posteriorcorrs$age[which.max(x)],4)
                     }))
age.max.corr <- t(age.max.corr)
age.max.corr.median <- median(age.max.corr) #median age #bayes
age.max.corr.CI <- quantile(age.max.corr, probs = c(0.025, 0.975)) #compute the credible interval based on all draws
age.max.corr.lower <- age.max.corr.CI[[1]]
age.max.corr.upper <- age.max.corr.CI[[2]]

age.max.corr.median

## [1] 15.0243

age.max.corr.CI

##    2.5%   97.5% 
## 14.6516 15.3224

age.max.corr <- as.data.frame(age.max.corr)
ggplot(age.max.corr, aes(x = V1)) + 
geom_histogram(bins = 12, fill = "#FFE5B1") + 
geom_vline(xintercept = 15.0243, color = "#ba275f", size = .5) +
labs(x="\nAge", y="Posterior Draw Number\n") +
theme_classic() +
theme(axis.text.x = element_text(size = 6, family = "Arial", color = c("black")), axis.text.y = element_blank(), axis.title = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
scale_x_continuous(breaks=c(14.5, 15, 15.5))   

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure4/AgeMaxCorr_histogram.pdf", device = "pdf", dpi = 500, width = 1.03 , height = 1.02)

Age of first null (zero) sensorimotor-association axis correlation: posterior median value + 95% credible interval

age.zero.corr <- deriv.SAaxis.posteriorcorrs %>% #find the age at which the correlation first hits 0 for each draw
    summarise(across(contains("draw"),
                     .fns = function(x){
                      round(min(deriv.SAaxis.posteriorcorrs$age[round(x, digits = 1) == 0]),4)
                     }))
age.zero.corr <- t(age.zero.corr)
age.zero.corr.median <- median(age.zero.corr) #median age #bayes
age.zero.corr.CI <- quantile(age.zero.corr, probs = c(0.025, 0.975)) #compute the credible interval based on all draws
age.zero.corr.lower <- age.zero.corr.CI[[1]]
age.zero.corr.upper <- age.zero.corr.CI[[2]]

age.zero.corr.median

## [1] 19.273

age.zero.corr.CI

##    2.5%   97.5% 
## 18.6767 20.1675

Sensorimotor-association axis correlation value at each age: posterior median values + 95% credible interval

deriv.SAaxis.posteriorcorrs <- deriv.SAaxis.posteriorcorrs %>% select(contains("draw"))

deriv.SAaxis.mediancorr <- apply(deriv.SAaxis.posteriorcorrs, 1, function(x){median(x)}) #median correlation value at each age

cor.credible.intervals <- apply(deriv.SAaxis.posteriorcorrs, 1, function(x){quantile(x, probs = c(0.025, 0.975))}) #compute the credible interval for the correlation value at each age based on all draws
cor.credible.intervals <- t(cor.credible.intervals)
cor.credible.intervals <- as.data.frame(cor.credible.intervals)
cor.credible.intervals <- cbind(cor.credible.intervals, (gam.derivatives.glasser %>% group_by(age) %>% group_keys()))
cor.credible.intervals$SA.correlation <- deriv.SAaxis.mediancorr 
colnames(cor.credible.intervals) <- c("lower","upper","age","SA.correlation")

Maximal correlation value

cor.credible.intervals %>% arrange(desc(SA.correlation)) %>% slice_head()

##       lower     upper      age SA.correlation
## 1 0.6629209 0.7040393 15.02429      0.6837556

Developmental alignment with the sensorimotor-association axis plot

cor.credible.intervals$max.corr.CI <- (cor.credible.intervals$age > age.max.corr.lower & cor.credible.intervals$age < age.max.corr.upper) #add a column that indicates whether each age is included in the age of maximal correlation credible interval (T/F)
cor.credible.intervals$max.cor.window <- cor.credible.intervals$age*cor.credible.intervals$max.corr.CI #add a column that only includes ages in this interval
cor.credible.intervals$max.cor.window[cor.credible.intervals$max.cor.window == 0] <- NA
                        
cor.credible.intervals$zero.corr.CI <- (cor.credible.intervals$lower < 0 & cor.credible.intervals$upper > 0) #add a column that indicates whether the credible interval for the correlation includes 0
cor.credible.intervals$zero.corr.window <- cor.credible.intervals$age*cor.credible.intervals$zero.corr.CI #add a column that only includes ages in the zero window
cor.credible.intervals$zero.corr.window[cor.credible.intervals$zero.corr.window == 0] <- NA

Figure 4B

ggplot(cor.credible.intervals, aes(x = age, y = SA.correlation, ymin = lower, ymax = upper)) + 
geom_line(size = .5) +
geom_ribbon(alpha = .2, fill = c("grey30")) +
labs(x="\nAge", y="Developmental Alignment to the Axis\n") +
geom_ribbon(aes(x = max.cor.window, y = SA.correlation), fill = "#ba275f") +
theme_classic() + 
theme(axis.text = element_text(size = 7, family = "Arial", color = c("black")), axis.title = element_text(size = 7, family = "Arial", color = c("black")), axis.line = element_line(size =.22), axis.ticks = element_line(size = .22)) +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .1)) +
scale_y_continuous(breaks=c(0,0.2,0.4,0.6)) 

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure4/Derivative_SAaxis_correlation_plot.pdf", device = "pdf", dpi = 500, width = 3.82 , height = 2.35)

Environmental influences on intrinsic fMRI activity vary across the sensorimotor-association axis during adolescence

environment <- read.csv("/cbica/projects/spatiotemp_dev_plasticity/sample_info/PNC/n1601_go1_environment_factor_scores_tymoore_20150909.csv")
environment <- environment %>% select(-bblid)
fluctuations.glasser.pnc <- merge(fluctuations.glasser.pnc, environment, by="scanid", sort = F)

Correlation between environment factor scores and age

cor.test(fluctuations.glasser.pnc$age, fluctuations.glasser.pnc$envSES)

## 
##  Pearson's product-moment correlation
## 
## data:  fluctuations.glasser.pnc$age and fluctuations.glasser.pnc$envSES
## t = 0.47644, df = 1031, p-value = 0.6339
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.04619975  0.07576231
## sample estimates:
##        cor 
## 0.01483646

Significant associations between fluctuation amplitude and environment factor scores

Number of cortical regions showing significant environment effects

sum(p.adjust(df.env$Anova.env.pvalue, method=c("fdr")) < 0.05) 

## [1] 141

• Number of cortical regions that still show a significant environment effect, after controlling for parental education

df.env %>% filter(p.adjust(df.env$Anova.env.pvalue, method=c("fdr")) < 0.05 & p.adjust(df.env$Anova.env.edu.pvalue, method=c("fdr")) < 0.05) %>% nrow()

## [1] 101

Percent of cortical regions showing significant environment effects

(sum(p.adjust(df.env$Anova.env.pvalue, method=c("fdr")) < 0.05))/(nrow(df.env))*100 

## [1] 41.96429

Across-cortex environment effect size and direction

ggseg(.data = df.env, atlas = "glasser", mapping=aes(fill=GAM.env.tvalue), position = c("stacked")) + 
  theme_void() + 
  labs(fill="Age Effect") +
  scale_fill_gradient2(low= "#F9AE76", mid = "white", high = "#952162", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0, limits=c(-5,5), oob=squish) +
  theme(legend.text = element_text(size = 6, family = "Arial", color = c("black")))

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure6/Brainmap_enveffect.jpg", device = "pdf", dpi = 500, width = 3.05 , height = 2.25)

Correlation between environment effects and the sensorimotor-association axis

All regions

cor.test(df.env$GAM.env.tvalue, df.env$SA.axis, method=c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  df.env$GAM.env.tvalue and df.env$SA.axis
## S = 3274572, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4820453

perm.sphere.p(df.env.spin$GAM.env.tvalue, df.env.spin$SA.axis, perm.id.full, corr.type='spearman')

## [1] 0

• All regions, controlling for parental education

cor.test(df.env$GAM.env.edu.tvalue, df.env$SA.axis, method=c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  df.env$GAM.env.edu.tvalue and df.env$SA.axis
## S = 2914840, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5389458

perm.sphere.p(df.env.spin$GAM.env.edu.tvalue, df.env.spin$SA.axis, perm.id.full, corr.type='spearman')

## [1] 0

Correlation plot

df.env$significant <- p.adjust(df.env$Anova.env.pvalue, method=c("fdr")) < 0.05

ggplot(df.env, aes(x = SA.axis, y = GAM.env.tvalue, fill = GAM.env.tvalue, color = significant)) + 
geom_point(size = 1.5, stroke = .12, pch = 21) +
scale_fill_gradient2(low= "#F9AE76", mid = "white", high = "#952162", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0, limits = c(-5,5), oob = squish) +
scale_color_manual(values = c("white","black")) +
labs(x="\nSensorimtor-Association Axis Rank", y="Environment Effect\n") +
geom_smooth(method = 'lm', se = TRUE, fill = alpha(c("gray70"),.7), col = "black", size = .5) +
theme_classic() + 
theme(legend.position = "none") +
theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size = .22), axis.ticks = element_line(size = .22)) +
scale_y_continuous(breaks = c(-4,  0, 4,  8))

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure6/AxisEnv_correlationplot.jpg", device = "pdf", dpi = 500, width = 2.5 , height = 1.7)

Predicted developmental trajectories for 10th and 90th percentile factor scores

SAaxismooths.by.env <- gam.fitted.env %>% group_by(envSES, age, SA.axis.bin) %>% do(est.mean = mean(.$fitted.centered)) %>% unnest(cols = c(est.mean)) #calculate the smooth fit for each S-A axis decile for low and high environment factor scores

SA.colorvector <- c("#FFC228", "#FFCA4E", "#FFD16A", "#FFDA89", "#FFE6B0", "#D7BCDA", "#BE93C4", "#9859A4", "#863E95", "#6F1382")

for(decile in c(1, 3, 6, 8, 10)){
  
  SAaxismooths.by.env.toplot <- SAaxismooths.by.env %>% filter(SA.axis.bin == decile)

  ggplot(data = SAaxismooths.by.env.toplot, aes(x = age, y = est.mean, group = interaction(envSES, SA.axis.bin), linetype = envSES )) +
  geom_line(size = .5, color = SA.colorvector[decile]) +
  theme_classic() +
  labs(x = "\nAge", y = "Predicted Fluctuation Amplitude\n") +
  theme(legend.position = "none")  +
  theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size=.22), axis.ticks = element_line(size = .22)) + scale_x_continuous(breaks=c(8, 12, 16, 20), expand = c(0,.45)) +
  scale_y_continuous(limits=c(-.06,.23), breaks=c(0,.1,.2))
  
  ggsave(filename = sprintf("/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure6/devtrajectories_env_decile%s.pdf", decile), device = "pdf", dpi = 500, width = 1.75 , height = 1.5) }

Environment effects in children, adolescent, and early adult samples

Correlations between environment t-values across developmental stages

cor.test(df.env.devstage$GAM.env.tvalue.child, df.env.devstage$GAM.env.tvalue.adolescent)

## 
##  Pearson's product-moment correlation
## 
## data:  df.env.devstage$GAM.env.tvalue.child and df.env.devstage$GAM.env.tvalue.adolescent
## t = 19.209, df = 334, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6693840 0.7716676
## sample estimates:
##       cor 
## 0.7244901

cor.test(df.env.devstage$GAM.env.tvalue.adolescent, df.env.devstage$GAM.env.tvalue.adult)

## 
##  Pearson's product-moment correlation
## 
## data:  df.env.devstage$GAM.env.tvalue.adolescent and df.env.devstage$GAM.env.tvalue.adult
## t = 17.989, df = 334, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6427370 0.7520399
## sample estimates:
##       cor 
## 0.7014904

cor.test(df.env.devstage$GAM.env.tvalue.child, df.env.devstage$GAM.env.tvalue.adult)

## 
##  Pearson's product-moment correlation
## 
## data:  df.env.devstage$GAM.env.tvalue.child and df.env.devstage$GAM.env.tvalue.adult
## t = 15.68, df = 334, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5849190 0.7087752
## sample estimates:
##       cor 
## 0.6511617

Correlation between environment effects and the sensorimotor-association axis in each developmental stage

Children

cor.test(df.env.devstage$GAM.env.tvalue.child, df.env.devstage$SA.axis, method=c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  df.env.devstage$GAM.env.tvalue.child and df.env.devstage$SA.axis
## S = 4323254, p-value = 3.778e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3161702

perm.sphere.p(df.env.devstage.spin$GAM.env.tvalue.child, df.env.devstage.spin$SA.axis, perm.id.full, corr.type='spearman')

## [1] 0.00765

Adolescents

cor.test(df.env.devstage$GAM.env.tvalue.adolescent, df.env.devstage$SA.axis, method=c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  df.env.devstage$GAM.env.tvalue.adolescent and df.env.devstage$SA.axis
## S = 2219600, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6489152

perm.sphere.p(df.env.devstage.spin$GAM.env.tvalue.adolescent, df.env.devstage.spin$SA.axis, perm.id.full, corr.type='spearman')

## [1] 0

Young Adults

cor.test(df.env.devstage$GAM.env.tvalue.adult, df.env.devstage$SA.axis, method=c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  df.env.devstage$GAM.env.tvalue.adult and df.env.devstage$SA.axis
## S = 4375798, p-value = 9.904e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.3078591

perm.sphere.p(df.env.devstage.spin$GAM.env.tvalue.adult, df.env.devstage.spin$SA.axis, perm.id.full, corr.type='spearman')

## [1] 0.0128

Across-cortex environment effect size and direction

Children

ggseg(.data = df.env.devstage, atlas = "glasser", mapping=aes(fill=GAM.env.tvalue.child), position = c("stacked")) + 
  theme_void() + 
  labs(fill="Age Effect") +
  scale_fill_gradient2(low= "#F9AE76", mid = "white", high = "#952162", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0, limits=c(-3,3), oob=squish) +
  theme(legend.text = element_text(size = 6, family = "Arial", color = c("black")))

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure6/Brainmap_enveffect_child.jpg", device = "pdf", dpi = 500, width = 2.5 , height = 2)

Adolescents

ggseg(.data = df.env.devstage, atlas = "glasser", mapping=aes(fill=GAM.env.tvalue.adolescent), position = c("stacked")) + 
  theme_void() + 
  labs(fill="Age Effect") +
  scale_fill_gradient2(low= "#F9AE76", mid = "white", high = "#952162", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0, limits=c(-3,3), oob=squish) +
  theme(legend.text = element_text(size = 6, family = "Arial", color = c("black")))

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure6/Brainmap_enveffect_adolescent.jpg", device = "pdf", dpi = 500, width = 2.5 , height = 2)

Young Adults

ggseg(.data = df.env.devstage, atlas = "glasser", mapping=aes(fill=GAM.env.tvalue.adult), position = c("stacked")) + 
  theme_void() + 
  labs(fill="Age Effect") +
  scale_fill_gradient2(low= "#F9AE76", mid = "white", high = "#952162", guide = "colourbar", aesthetics = "fill", name = NULL, midpoint = 0, limits=c(-3,3), oob=squish) +
  theme(legend.text = element_text(size = 6, family = "Arial", color = c("black")))

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure6/Brainmap_enveffect_adult.jpg", device = "pdf", dpi = 500, width = 2.5 , height = 2)

Age-specific correlations between environment effects and the S-A Axis

env.SAaxis.posteriorcorrs <- read.table("/cbica/projects/spatiotemp_dev_plasticity/FluctuationAmplitude/GAMRESULTS/SAaxis_environmenteffects_correlation_byage_glasser.csv", header = F, sep = ",") #get the regional derivative-regional axis rank correlation at each age for all draws from the posterior distribution
colnames(env.SAaxis.posteriorcorrs) <- sprintf("draw%s",seq(from = 1, to = ncol(env.SAaxis.posteriorcorrs)))
env.SAaxis.posteriorcorrs$age <- seq(min(fluctuations.glasser.pnc$age), max(fluctuations.glasser.pnc$age), length.out = 200)

Age of maximal sensorimotor-association axis correlation: posterior median value + 95% credible interval

age.max.corr <- env.SAaxis.posteriorcorrs %>% #find the age at which the correlation is largest for each draw
    summarise(across(contains("draw"),
                     .fns = function(x){
                       round(env.SAaxis.posteriorcorrs$age[which.max(x)],4)
                     }))
age.max.corr <- t(age.max.corr)
age.max.corr.median <- median(age.max.corr) #median age #bayes
age.max.corr.CI <- quantile(age.max.corr, probs = c(0.025, 0.975)) #compute the credible interval based on all draws
age.max.corr.lower <- age.max.corr.CI[[1]]
age.max.corr.upper <- age.max.corr.CI[[2]]

age.max.corr.median

## [1] 15.5461

age.max.corr.CI

##    2.5%   97.5% 
## 14.4280 16.5896

age.max.corr <- as.data.frame(age.max.corr)
ggplot(age.max.corr, aes(x = V1)) + 
geom_histogram(bins = 10, fill = "#FFE5B1") + 
geom_vline(xintercept =  age.max.corr.median, color = "#ba275f", size = .5) +
labs(x="\nAge", y="Posterior Draw Number\n") +
theme_classic() +
theme(axis.text.x = element_text(size = 6, family = "Arial", color = c("black")), axis.text.y = element_blank(), axis.title = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
scale_x_continuous(breaks=c(14, 15, 16, 17))   

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure6/AgeMaxCorr_histogram.pdf", device = "pdf", dpi = 500, width = .7 , height = .7)

Sensorimotor-association axis correlation value at each age: posterior median values + 95% credible interval

env.SAaxis.posteriorcorrs <- env.SAaxis.posteriorcorrs %>% select(contains("draw"))

env.SAaxis.mediancorr <- apply(env.SAaxis.posteriorcorrs, 1, function(x){median(x)}) #median correlation value at each age

cor.credible.intervals.env <- apply(env.SAaxis.posteriorcorrs, 1, function(x){quantile(x, probs = c(0.025, 0.975))}) #compute the credible interval for the correlation value at each age based on all draws
cor.credible.intervals.env <- t(cor.credible.intervals.env)
cor.credible.intervals.env <- as.data.frame(cor.credible.intervals.env)
cor.credible.intervals.env$age <- seq(min(fluctuations.glasser.pnc$age), max(fluctuations.glasser.pnc$age), length.out = 200)
cor.credible.intervals.env$SA.correlation <- env.SAaxis.mediancorr 
colnames(cor.credible.intervals.env) <- c("lower","upper","age","SA.correlation")

Maximal correlation value

cor.credible.intervals.env %>% arrange(desc(SA.correlation)) %>% slice_head()

##       lower     upper     age SA.correlation
## 1 0.4609763 0.5500586 15.6206      0.5064657

Environment effects along the sensorimotor-association axis plot

cor.credible.intervals.env$max.corr.CI <- (cor.credible.intervals.env$age > age.max.corr.lower & cor.credible.intervals.env$age < age.max.corr.upper) #add a column that indicates whether each age is included in the age of maximal correlation credible interval (T/F)
cor.credible.intervals.env$max.cor.window <- cor.credible.intervals.env$age*cor.credible.intervals.env$max.corr.CI #add a column that only includes ages in this interval
cor.credible.intervals.env$max.cor.window[cor.credible.intervals.env$max.cor.window == 0] <- NA

cor.credible.intervals.env$zero.corr.CI <- (cor.credible.intervals.env$lower < 0 & cor.credible.intervals.env$upper > 0) #add a column that indicates whether the credible interval for the correlation includes 0
cor.credible.intervals.env$zero.corr.window <- cor.credible.intervals.env$age*cor.credible.intervals.env$zero.corr.CI #add a column that only includes ages in the zero window
cor.credible.intervals.env$zero.corr.window[cor.credible.intervals.env$zero.corr.window == 0] <- NA
cor.credible.intervals.env$zero.corr.window.young <- NA
cor.credible.intervals.env$zero.corr.window.young[1:7] <- cor.credible.intervals.env$zero.corr.window[1:7]
cor.credible.intervals.env$zero.corr.window.old <- NA
cor.credible.intervals.env$zero.corr.window.old[174:200] <- cor.credible.intervals.env$zero.corr.window[174:200]

ggplot(cor.credible.intervals.env, aes(x = age, y = SA.correlation, ymin = lower, ymax = upper)) + 
geom_line(size = .4) +
geom_ribbon(alpha = .2, fill = c("grey30")) +
labs(x="\nAge", y="Environment Effect Correlation with the S-A Axis\n") +
geom_ribbon(aes(x = max.cor.window, y = SA.correlation), fill = "#F9AE76") +
geom_ribbon(aes(x = zero.corr.window.young, y = SA.correlation), fill = "gray60") +
geom_ribbon(aes(x = zero.corr.window.old, y = SA.correlation), fill = "gray60") +
theme_classic() + 
theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size =.22), axis.ticks = element_line(size = .22)) +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .1)) +
scale_y_continuous(breaks=c(-.2,0,0.2,0.4,0.6), limits=c(-.2,0.6)) 

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure6/EnvEffect_SAaxis_age_correlations.jpg", device = "pdf", dpi = 500, width = 2.4 , height = 1.6)

Region- and age-specific environment effect magnitude

Region 3b

rh_R_3b.fit <- gam.varyingcoefficients(measure = "fluctuations", atlas = "glasser", dataset = "pnc",region = "rh_R_3b", smooth_var = "age", int_var = "envSES", covariates = "sex + RMSmotion", knots = 3, set_fx = TRUE, increments = 200, draws = 1, return_posterior_coefficients = FALSE)

Region 9m

rh_R_9m.fit <- gam.varyingcoefficients(measure = "fluctuations", atlas = "glasser", dataset = "pnc",region = "rh_R_9m", smooth_var = "age", int_var = "envSES", covariates = "sex + RMSmotion", knots = 3, set_fx = TRUE, increments = 200, draws = 1, return_posterior_coefficients = FALSE)

ggplot() +
geom_line(data = rh_R_3b.fit, aes(x = age, y = envSES.slope), size = .75, color = c("#FFD16A")) +
geom_line(data = rh_R_9m.fit, aes(x = age, y = envSES.slope), size = .75, color = c("#9859A4")) +
labs(x="\nAge", y="Regional Environment Effect\n") +
theme_classic() + 
theme(axis.text = element_text(size = 6, family = "Arial", color = c("black")), axis.title = element_text(size = 6, family = "Arial", color = c("black")), axis.line = element_line(size =.22), axis.ticks = element_line(size = .22)) +
scale_x_continuous(breaks=c(8, 12, 16, 20), expand = c(0,.45))  +
scale_y_continuous(breaks=c(-.02,0,0.02)) +
geom_hline(yintercept = 0, linetype = "dotted", color = "grey70", size = .4)

ggsave(filename = "/cbica/projects/spatiotemp_dev_plasticity/figures/Manuscript_Figures/NatNeuro/Figure6/Regional_ageresolved_effects.jpg", device = "pdf", dpi = 500, width = 1.9 , height = 1.6)